[WIP] Restraints API #1043
[WIP] Restraints API #1043
Conversation
|
🚨 API breaking changes detected! 🚨 |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1043 +/- ##
==========================================
- Coverage 94.50% 90.03% -4.47%
==========================================
Files 134 147 +13
Lines 9986 10689 +703
==========================================
+ Hits 9437 9624 +187
- Misses 549 1065 +516
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
|
🚨 API breaking changes detected! 🚨 |
| ] | ||
| for at3 in atom_pool: | ||
| if at3 != atom_idx and at3 in at2_neighbors: | ||
| angles.append((atom_idx, at2, at3)) |
There was a problem hiding this comment.
From today's discussion: we need to avoid cases where we can cross rings with a non-aromatic bond. Most likely easiest way to do this is to have a special case here that can optionally filter by ring (i.e. only pick angles where you have atoms indices that exist solely within one ring).
| raise ValueError(errmsg) | ||
|
|
||
| # 2. Get the central atom | ||
| center = get_central_atom_idx(rdmol) |
There was a problem hiding this comment.
We may need to use the old approach - need to test this and see if nx is picking long paths when it comes to rings.
| bond, ang1, ang2, dih1, dih2, dih3 = _get_restraint_distances( | ||
| atomgroup | ||
| ) | ||
|
|
There was a problem hiding this comment.
TODO: add in checks here so we can warn if it's a bad pick. May push this to a follow-up issue.
| topology_format=_get_mda_topology_format(topology), | ||
| ) | ||
|
|
||
| host_ag1 = u.atoms[host_idxs] |
There was a problem hiding this comment.
Maybe find better words for this (or maybe actually explain what the variable means!)
| Optional[list[int]] | ||
| A list of indices for a selected combination of H0, H1, and H2. | ||
| """ | ||
| h0_eval = EvaluateHostAtoms1( |
There was a problem hiding this comment.
Maybe add some comments on what is going on here?
|
|
||
| Note | ||
| ---- | ||
| We assume the temperature to be 298.15 Kelvin. |
There was a problem hiding this comment.
Add acknowledgements to Vytas.
| for i, at in enumerate(self.host_atom_pool): | ||
| distance_bounds = all(self.results.distances[i] > self.minimum_distance) | ||
| mean_angle = circmean(self.results.angles[i], high=np.pi, low=0) | ||
| angle_bounds = check_angle_not_flat( |
There was a problem hiding this comment.
Let's try doing this for every frame.
| width=1.745 * unit.radians, | ||
| ) | ||
| mean_dihed = circmean(self.results.dihedrals[i], high=np.pi, low=-np.pi) | ||
| dihed_bounds = check_dihedral_bounds(mean_dihed) |
| not_collinear, | ||
| ] | ||
| ): | ||
| self.results.valid[i] = True |
There was a problem hiding this comment.
If you do that, then you need to set things to False too!
| for i, valid_h0 in enumerate(h0_eval.results.valid): | ||
| if valid_h0: | ||
| g1g2h0_atoms = guest_atoms.atoms[1:] + host_atom_pool.atoms[i] | ||
| h1_eval = EvaluateHostAtoms1( |
There was a problem hiding this comment.
TODO: fix the order, this is incorrect.
| h1_eval = EvaluateHostAtoms1( | ||
| g1g2h0_atoms, | ||
| host_atom_pool, | ||
| minimum_distance, |
There was a problem hiding this comment.
Let's double check this is indeed 0.5 nm (or higher)
| temperature, | ||
| ) | ||
| for j, valid_h1 in enumerate(h1_eval.results.valid): | ||
| g2h0h1_atoms = g1g2h0_atoms.atoms[1:] + host_atom_pool.atoms[j] |
| distance2_bounds = all(self.results.distances2[i] > self.minimum_distance) | ||
| mean_dihed = circmean(self.results.dihedrals[i], high=np.pi, low=-np.pi) | ||
| dihed_bounds = check_dihedral_bounds(mean_dihed) | ||
| dihed_variance = check_angular_variance( |
There was a problem hiding this comment.
Same here could do per frame.
| not_collinear, | ||
| ] | ||
| ): | ||
| self.results.valid[i] = True |
There was a problem hiding this comment.
Here again, make sure you fill in the False case.
| ) | ||
|
|
||
| if any(h2_eval.ressults.valid): | ||
| d1_avgs = np.array([d.mean() for d in h2_eval.results.distances1]) |
There was a problem hiding this comment.
This needs to be filtered by True!
| raise ValueError(errmsg) | ||
|
|
||
| # Set the equilibrium values as those of the final frame | ||
| u.trajectory[-1] |
There was a problem hiding this comment.
Raise an issue / todo here to experiment on what frame / values to pick - ideally some kind of mean with the frame closest to the mean values.
There was a problem hiding this comment.
Thanks @IAlibay , just a few small comments from a first quick pass. I left out the boresch.py and the utils script for now as discussed. Will do another pass after the break!
| Partially reproduced from Yank. | ||
| """ | ||
|
|
||
| lambda_restraints = GlobalParameterState.GlobalParameter( |
There was a problem hiding this comment.
Does 1 here mean that the restraint is fully turned on? Maybe add a comment here?
There was a problem hiding this comment.
Took me a while to remember, the answer is yes - assuming the underlying Force does it properly ;)
What I mean here is that technically it's up to the underlying for to define how the controlling parameter is applied, so in all the cases below it's "lambda_restraints * E" but nothing really stops it from being "1/lambda_restraints * E" is someone wanted to cause chaos.
That being said, I'm pretty sure we'd do our due dilligence to not let that happen..
This is all to say "Yes I've added the comment and updated the docstring" :)
|
No API break detected ✅ |
Checklist
newsentryDevelopers certificate of origin